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Abstract By varying the noise intensity, we study stochastic spiking coherence (i.e., col- 
lective coherence between noise-induced neural spikings) in an inhibitory population of 
subthreshold neurons (which cannot fire spontaneously without noise). This stochastic spik- 
ing coherence may be well visualized in the raster plot of neural spikes. For a coherent case, 
partially-occupied "stripes" (composed of spikes and indicating collective coherence) are 
formed in the raster plot. This partial occupation occurs due to "stochastic spike skipping" 
which is well shown in the multi-peaked interspike interval histogram. The main purpose 
of our work is to quantitatively measure the degree of stochastic spiking coherence seen in 
the raster plot. We introduce a new spike-based coherence measure by considering the 
occupation pattern and the pacing pattern of spikes in the stripes. In particular, the pacing 
degree between spikes is determined in a statistical-mechanical way by quantifying the aver- 
age contribution of (microscopic) individual spikes to the (macroscopic) ensemble-averaged 
global potential. This "statistical-mechanical" measure M, is in contrast to the conventional 
measures such as the "thermodynamic" order parameter (which concerns the time-averaged 
fluctuations of the macroscopic global potential), the "microscopic" correlation-based mea- 
sure (based on the cross-correlation between the microscopic individual potentials), and the 
measures of precise spike timing (based on the peri-stimulus time histogram). In terms of 
Ms, we quantitatively characterize the stochastic spiking coherence, and find that M, re- 
flects the degree of collective spiking coherence seen in the raster plot very well. Hence, 
the "statistical-mechanical" spike-based measure may be used usefully to quantify the 
degree of stochastic spiking coherence in a statistical-mechanical way. 

Keywords Inhibitory Subthreshold Neurons ■ Stochastic Spiking Coherence ■ Statistical- 
Mechanical Measure 



1 Introduction 

Recently, brain rhythms have attracted much attention teuzsaki 1 1200^ . Synchronous oscil- 
lations in neural systems may be used for efficient sensory processing (e.g., binding of the 
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integrat ed whole ima ge in the visual cortex is accomplished via synchronization of neural 
firings) ( lGravlll994h . In addition to such neural encoding of sensory stimuli, neural syn- 
chronization is also correlated with pathological rhythms associa ted with neural diseases 
(e.g., epileptic seizures and tremors in the Parkinson's disease) l lHauschildtetal.ll2006l ; 
iGrosse et al. Il2003 : iTass et ai7lll998h . Here, we are interested in these synchronized neu- 
ral oscillations. 

A neural circuit in the major parts of the brain such as thalamus, hippocampus, and 
cortex consists of a few types of excitatory principal cells and diverse types of inhibitory in- 
temeurons. Functional diversity of intemeurons i ncreases the computational power of prin- 
cipal cells ( iBuzsaki 1 12OO6I : iBuzsaki et al. II2004I) . To understand the mechanisms of syn- 
chronous brain rhythms, neural systems composed of excitatory neurons and/or inhibitory 
neurons hav e been investi gated, and thus three types of synchronization mechanisms have 
been found (IWang||2003h : recurrent excitation between principal cells, mutual inhibition 
between interneurons, and feedback between excitatory and inhibitory ne urons. Perfect syn- 
chronization occurs in the network of pulse-coupled excitatory neurons jKuramotoll 199lL 
iMiroUo & Strogatz1ll990h . However, this synchronization via mutual excitation cannot be 
always stable in networks with slowly decaying synaptic couplings. Stability of synchro- 
nization is shown to depend on both the time course of synaptic interaction and t he response 
of neurons to small depolarization jVreeswiik et al. I [l99?: 'Hanse l et al. I [l993). When the 
decay time of the synaptic interaction is enough long, mutual inhibition (rather than exci- 
tation) may synchronize neural firing activities. By providing a coherent oscillatory output 
to the principal cells, interneuronal networks play the role of the backbones (i.e., synchro- 
nizers or pac emakers) of many bra in oscillations such as the 10-Hz thalamocortical spin- 
dle rhythms dWang & Rinzel Ill992l : lOolomb & Rinzellll994) and the 40-Hz fast gamma 



aie rnytnms a Wang & KmzeiiiiWii : luoiomb & Kmzeiwm^) and tne 4U-Hz tast gamma 
rhythms in the hippocampus and the neocortex dWang & BuzsaM 1ll99d : IWhiteetal.il 19981 : 
IWhittington et al. ibood : iTiesinga et al. 1 l200lh . When the feedback between the excitatory 
and inhibitory popul ations is strong, neural synchrony appears via the "cross-talk" between 
the two populations Iwhittington et al. 1120001 : iTiesinga et ^1120011 : iHansel & Mato II2003I : 
iBorgers &Kopellll2003l . 120051) . 



Most past studies exploring mechanisms of neural synchronization were done in neural 
systems composed of spontaneously firing (i.e., self-oscillating) suprathreshold neurons. For 
this case, neural coherence occurs via cooperation of regular firings of suprathreshold self- 
firing neurons. In contrast, neural systems composed of subthreshold neurons have received 
little attention. Unlike the suprathreshold case, each subthreshold neuron in the absence of 
coupling cannot fire spontaneously without noise; it can fire only with the help of noise. 
Recently, stochastic spiking coherence (i.e., collective coherence emerging via cooperation 
of noise-induced s pikings) was observed in an exc itatory population of pulse-coupled sub- 
threshold neurons (IWang et al.]|2000l : lLim & Kimll2 007. 2009). This kind of works may be 
thought to correspond to a "subthreshold version" of neural synchronization through mutual 
excitation. Due to the stochastic spiking coherence, synaptic current, injected into each in- 
dividual neuron, becomes temporally coherent. Hence, temporal coherence resonance of an 
individual subthreshold neuron in the network may be enhanced. This enhancement of co- 
herence resonance in an excitatory network of subthreshold Hodgkin-Huxley neurons was 
characterized in terms of the coherence factor jS, repre senting the degree of sharpness of 
the peak in the power spectrum of an individual neuron ('Wangetal."200d). In this way, the 
measure /3 is used for characterization of the temporal coherence of an individual neuron. 
However, we note that jS is not a measure for directly measuring the degree of collective 
coherence in the whole population. 
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In this paper, we are concerned about the "subthreshold version" of neural synchro- 
nization via mutual inhibition. In Section |2l we describe the biolog ical conductance-based 
Morris-Lecar (ML) neuron model with voltage-ga ted ion channels jMorris & LecarlflQEll ; 
iRinzel & Ermentrout 1 1 19981 : iTsumoto et al. |[200a) . The ML neurons (used in our study) 
exhibit the type-II excitability (i.e., the firing frequency begins to increase from a non- 
zero value when the stimulus exceeds a threshold value), and they interact via inhibitory 
GABAergic synapses whose activity increases fast and decays slowly. In Section [3] we 
characterize stochastic spiking coherence in a large population of inhibitory subthreshold 
ML neurons by varying the noise intensity for a fixed coupling strength. Weakly coherent 
states with oscillating ensemble-averaged global membrane potential Vq are thus found to 
appear in a range of intermediate noise intensity [i.e., regular global oscillation (with re- 
duced amplitude and increased frequency) emerges via cooperation of irregular individual 
oscillations]. Emergence of collective coherence may be well described in terms of the con- 
ventional "thermodynamic" order parameter which concerns the time-averaged fluctuations 
of the macroscopic global potential Vc- We note that this stochastic spiking coherence may 
be well visualized in the raster plot of neural spikes (i.e., a spatiotemporal plot of neural 
spikes) which is directly obtained in experiments. For the coherent case, "stripes" (com- 
posed of spikes and indicating collective coherence) are found to be formed in the raster 
plot. Due to coherent contribution of spikes, local maxima of the global potential Vc appear 
at the centers of stripes. However, these stripes are partially occupied. Individual inhibitory 
neurons exhibit intermittent spikings phase-locked to Vc at random multiples of the period of 
Vc- This "stochastic phase locking" leading to "stochastic spike skipping" is well shown in 
the interspike interval (ISI) histogram with multiple peaks. The multi-peaked ISI histogram 
shows some indication of weak collective spiking coherence. 

The main purpose of our work is to quantitatively measure the degree of stochastic spik- 
ing coherence seen in the raster plot. In Section [3l we introduce a new type of spike-based 
coherence measure by taking into consideration the occupation pattern and the pacing 
pattern of spikes in the stripes of the raster plot. In particular, the pacing degree between 
spikes is determined in a statistical-mechanical way by quantifying the average contribution 
of (microscopic) individual spikes to the (macroscopic) global potential Vc - This "statistical- 
mechanical" measure is in contrast to the conven t ional measures such as the "thermo- 
dynamic" order parameter l lGolomb & Rinzel" ] 1 1994 |Hansel&Mato]|2003|), the "micro- 
scopic" correlation-ba sed measure (based on the cross-correlation s between the microscopic 
individual potentials) dWang & Buzsaki ll 19961 : 1 White et al. Ill998h . a nd the measures of pre- 
cise spike timing based on the peri-stimulus time histogram (PSTH) llMainen & Seinowski I 
0*995; Schreiber et al. 2003). The "thermodynamic" order parameter and the "microscopic" 
measure concern just the the macroscopic global potential Vc and the microscopic individual 
potentials, respectively without considering any quantitative relation between Vc and the mi- 
croscopic individual potentials. (The auto-correlation of the global activity used in the work 
of Brunei and Hakim (1999) may also be regarded as a kind of "thermodynamic" measure.) 
For the PSTH-based measure "events," corresponding to peaks of the instantaneous popula- 
tion firing rate, are selected through setting a threshold. Then, the measures for the reliability 
and the precision of spike timing concern only the spikes within the events, in contrast to the 
case of the "statistical-mechanical" measure where all spikes are considered (without select- 
ing events). A main difference between the conventional and the new spike-based measures 
lies in determining the pacing degree of spikes. The precision of spike timing for the con- 
ventional case is given by just the standard deviation of (microscopic) individual spike times 
within an event without considering the quantitative contribution of (microscopic) individ- 
ual spikes to the (macroscopic) global activity. Hence, the PSTH-based measure is not a 
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statistical-mechanical measure. However, if we take the instantaneous population firing rate 
as a global activity and exactly define the global cycles [see Fig.|4ja)] and the global phases 
[see Eqs. ( 115b and 1 II6H like our case, then the conventional PSTH-based measure may 
also develop into a similar statistical-mechanical measure. By varying the noise intensity, 
we quantitatively characterize the stochastic spiking coherence in terms of the "statistical- 
mechanical" measure Ms, and find that Mj reflects the degree of collective spiking coherence 
seen in the raster plot very well. We also expect that may be implemented for character- 
izing the degree of collective coherence in the experimentally-obtained raster plot of neural 
spikes. Finally, a summary is given in Section]?] 



2 Morris-Lecar Neuron Model 

In this section, we describe the biological neuron model used in our computational study. 
We consider an ensemble of globally coupled neurons. As an element in our neural sys- 
tem, we choose the conductance-based ML neuron model with voltage-gated ion channels, 
originally proposed to describe the ti me-evolution pattern o f the membrane potential for 
the giant m uscle fibers of barnacles ( iMorris & Lecar I ll98lL iRinzel & Ermentrout 1 1 19981 : 
IXsumoto et a l. 2006). The population dynamics in this neural network is governed by a set 
of the following differential equations: 

dvi 

C— = —Ihm.i + IdC + D^i — Isyn,i , (1) 



dt 

dwj _ {w„{vi)-Wi) 

dt Tr(v,) 

dsi 
dt 



(2) 

^ = a.s„{vi){l-si)-p.si, i=l,---,N, (3) 



where 



hon.i = ka.i+lK.i+h.i (4) 
= gCfl'«~(vi)(v; - Vca)+gKWi{Vi - VK)+gL{Vi - Vi), (5) 

hyn,i = J^—^ ^j{t){vi-V,,yn), (6) 

moo(v) = 0.5[l+tanh{(v-Vi)/y2}], (7) 
w„(v) = 0.5[l+tanh{(v-V3)/y4}], (8) 

Tr(v) = l/cosh{(v-V3)/(2V4)}, (9) 
s^{vi) = l/[l+e-(^'-"*'/^]. (10) 

Here, the state of the rth neuron at a time f (measured in units of ms) is characterized by 
three state variables: the membrane potential v, (measured in units of mV), the slow recov- 
ery variable w, representing the activation of the current (i.e., the fraction of open 
channels), and the synaptic gate variable Sj denoting the fraction of open synaptic ion chan- 
nels. In Eq. (O, C represents the capacitance of the membrane of each neuron, and the time 
evolution of v, is governed by four kinds of source currents. 

The total ionic current 7,o„., of the /'th neuron consists of the calcium current Ica.i, the 
potassium current /jf,;, and the leakage current //,.,. Each ionic current obeys Ohm's law. The 
constants gca, gK, and gi are the maximum conductances for the ion and leakage channels. 
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and the constants Vca, Vk, and Vi are the reversal potentials at which each current is bal- 
anced by the ionic concentration difference across the membrane. Since the calcium current 
Ica,i changes much faster than the potassium current Ikj, the gate variable m, for the Ca^+ 
channel is assumed to always take its saturation value m„{vi). On the other hand, the activa- 
tion variable w, for the channel approaches its saturation value Woo(v,) with a relaxation 
time TR{vi)/(j), where Tr has a dimension of ms and (j) is a (dimensionless) temperature-like 
time scale factor. 



Each ML neuron is also stimulated by the common DC current loc and an indepen- 
dent Gaussian white noise [see the 2nd and 3rd terms in Eq. lO] satisfying (<§, (r)) = and 
(^i{t) ^j{t')) = 5ij 5{t — t'), where (• • ■) denotes the ensemble average. The noise ^ is a para- 
metric one which randomly perturbs the strength of the applied current loc, ™d its intensity 
is controlled by the parameter D. The last term in Eq. l[T]l represents the coupling of the 
network. Each neuron is connected to all the other ones through global synaptic couplings. 
/„.„,, of Eq. ^ represents such synaptic current injected into the ith neuron. Here the cou- 
pling strength is controlled by the parameter / and Vjvn is the synaptic reversal potential. We 
use Vsxn = — 80 mV for the inhibitory synapse and Vjvn = mV for the excitatory synapse . 
The synaptic gate variab le s obeys the 1st order kinetics of Eq. (O ( Golomb & Rinzel 1 19941 : 
IWang&Buzsakilll996h . Here, the normalized concentration of synaptic transmitters, acti- 
vating the synapse, is assumed to be an instantaneous sigmoidal function of the membrane 
potential with a threshold v* in Eq. dlOb . where we set v* = mV and 5 = 2 mV. The 
transmitter release occurs only when the neuron emits a spike (i.e., its potential v is larger 
than V*). The synaptic channel opening rate, corresponding to the inverse of the synaptic rise 
time T,-, is a = 10 ms^^ which assures afastrise of /^v,, (Borge rs & Kope ll 2003, 2005). On 
the other hand, the synaptic closing rate jS , which is the inverse of the synaptic decay time 
Tij, depends on the type of the synaptic receptors. For the inhibitory GABAergic synapse 
(involving the GABAa receptors) we set jS = 0.1 ms^', and for the excitatory glutamate 
synapse (involving the AMPA receptors), we use jS = 0.5 ms" ( iBorgers & Kopell Il2003l 
|2005|). Thus, Isyn decays slowly. 



The ML ne uron may exhibit either type- I or type-II excitability, depending on the sys- 
tem parameters jRinzel & Ermentrout Ill998h . For the case of type-I (type-II), the firing fre- 
quency begins to increase from zero (a non-zero finite value) when loc passes a threshold. 
Throughout this paper, we consider the case of type-II excitability where gca =4.4 mS /cm? , gK = 
8 mS/cm^, = 2 mS/cm^, Vca = 120 mV, = -84 mV, Vl = -60 mV, C = 20 /iF/cm^, 
^ = 0.04, Vl = - 1.2 mV, ^2 = 18 mV, V3 = 2 mV, and Va = 30 mV. As loc passes a thresh- 
old in the absence of noise, each single type-II ML neuron begins to fire with a nonzero fre- 
quency that is relatively insensitive to the change in loc (Hodgkin 1948; Izhikevich 20 001). 
Numerical integration of Eqs. l[T][3]l is done using the Heun method ( San Miguel & Toral I 
I2OOOI) (with the time step At = 0.01 ms) similar to the second-order Runge-Kutta method, 
and data for (v,, w;,j,) (i = 1, . . . ,A^) are obtained with the sampling time interval At = \ ms. 
For each realization of the stochastic process in Eqs. l[T][3]l, we choose a random initial point 
[v,(0),w,(0),j,(0)] for the ith (/ = 1, . . . ,N) neuron with uniform probability in the range of 
v,(0) e (-70,50), w,(0) e (0.0,0.6), andj,(0) e (0.0,1.0). 
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Fig. 1 Intermittent noise-induced firings of a single subthreshold ML neuron: (a) time series of the membrane 
potential v and (b) interspike interval (ISI) histogram for Idc = 87 /jA/cm^ and Z) = 20 • ms'/^/ cm ; the 
ISI histogram is made of 5 x lO'' ISIs and the bin size for the histogram is 5 ms. Regular firings of a single 
suprathreshold ML neuron: (c) time series of v for loc = 95 /iA/cm^ and Z) = 0. 



3 Characterization of Stochastic Spiking Coherence in An Inhibitory Population of 
Subthreshold ML Neurons 

In this section, we study stochastic spiking coherence in a population of inhibitory sub- 
threshold ML neurons. The stochastic spiking coherence is quantitatively characterized in 
terms of a new "statistical-mechanical" spike-based measure. 

We first consider the case of a single ML neuron. An isolated subthreshold neuron cannot 
fire spontaneously without noise; it may generate firings only with the aid of noise. Figure 
[TJa) shows a time series of the membrane potential v of a subthreshold neuron for Ijx: = 87 
jUA/cm^ and D = 20 jUA ■ ms'/^/cm^. Noise-induced spikings appear intermittently. For 
this subthreshold case, the ISI histogram is shown in Fig. \V[b). The most probable value of 
the ISIs (corresponding to the main highest peak) is 97.5 ms. But, due to a long tail in the 
ISI histogram the average value of the ISIs becomes 161.6 ms. This noise-induced irregular 
oscillation of v is in contrast to the regular oscillation of v for a suprathreshold neuron (which 
fires spontaneously without noise) when loc = 95 fi A/ cm^ and D = Q [see Fig. [He)]. 

We consider an inhibitory population of A' globally coupled subthreshold ML neurons 
for loc = 87 jUA/cm^ and set the coupling strength as / = 3 mS/cm-. (Hereafter, for con- 
venience we omit the dimensions of loc, D, and J.) By varying the noise intensity D, we 
investigate the stochastic spiking coherence. Emergence of collective spiking coherence may 
be well described by the (population-averaged) global potential, 

Vcit) = ^tviit). (11) 

In the thermodynamic limit (A^ — > o°), a collective state becomes coherent if AVc{t) (= 
^c(0 ~ ^cit)) is non-stationary (i.e., an oscillating global potential Vc appears for a co- 
herent case), where the overbar represents the time average. Otherwise (i.e., when AVq is 
stationary), it becomes incoherent. Thus, the mean square deviation of the global potential 
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Fig. 2 Order parameter and partial occupation in N{= lO') globally coupled inhibitory subthreshold ML 
neurons: (a) plots of the order parameter versus the noise intensity D. (b) partially-occupied raster plot 
of spikings a ((': neuron index) and plot of number of spikes (N,) versus time (t) and (c) time series of the 
ensemble-averaged global potential Vq. Full occupation in N{= 10^ ) globally coupled excitatory ML neurons: 
(d) fully-occupied raster plot of spikings and plot of N, versus I and (e) time series of Va- For both inhibitory 
and excitatory cases, loc = 87 ^A/cm^, D = 20 /JA • ms'''^/cm^, 7=3 mS/cm^, and the bin size for in 
(b) and (d) is 1 ms 



Vc (i.e., time-averaged fluctuations of Vg), 



ff={VG{t)-VG{t))^ (12) 

plays t he role of an order pa rameter used for describing the coherence-incoherence tran- 
sition ( iManrubia et al. Il2004h . For the coherent (incoherent) state, the order parameter 
approaches a non-zero (zero) limit value as A' goes to the infinity. Figure [2 a) shows a plot 
of the order parameter versus the noise intensity. For D < D*i {~ 9.4), incoherent states exist 
because the order parameter & tends to zero as A' — > As D passes the lower threshold Djf, 
a coherent transition occurs because of a constructive role of noise to stimulate coherence 
between noise-induced spikings. However, for large D > DJ] (~ 33.4) such coherent states 
disappear (i.e., a transition to an incoherent state occurs when D passes the higher threshold 
Dp due to a destructive role of noise to spoil the collective spiking coherence. 
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As an example, we consider a coherent state for D = 20. For this case stochastic spiking 
coherence may be well visualized in terms of the raster plot of neural spikes. The upper 
panel in Fig.|2lb) shows the raster plot for A' = lO-' . We note that stripes (composed of spikes 
and indicating collective coherence) appear successively in the raster plot. The number of 
spikes (Ns) in each stripe is given in the lower panel. About 100 spikes constitute a stripe 
[i.e., only a fraction (about 1/10) of the total neurons fire in each stripe]. In this way partial 
occupation occurs in the stripes. A regularly oscillating global potential Vc emerges through 
cooperation of spikes in partially-occupied stripes. A time series of Vc is shown in Fig.|2tc). 
The amplitude of Vc is much smaller than that in Fig. [TJa) for the isolated single case. 
This reduction of amplitude occurs mainly because of partial occupation. Local maxima 
of Vc appear at the centers of stripes where the spiking densities are locally highest. The 
global period Tc of Vc (corresponding to the average intermax interval of Vc or equivalently 
corresponding to the average interstripe interval in the raster plot) is 54.2 ms. Hence the 
global frequency fa (= 18.5 Hz) of Vc is roughly as twice as the most probable frequency 
(= 10.3 Hz) of V for the isolated single case. Thus, a regular global oscillation with reduced 
amplitude and increased frequency occurs for the partially-occupied case. For comparison, 
we also consider the case of full occupation which occurs in a population of excitatory 
neurons for the same set of parameters foe, D, and /. [We emphasize that partial occupation 
may also occur even for the excitatory case of smaller J (e.g., / = 1).] Figures [2ld) and 
|2e) show the raster plot and a time series of Vc for the fully-occupied case, respectively. 
Fully-occupied stripes appear successively at nearly regular time intervals At (= 97.9 ms), 
in contrast to the partially-occupied case. A regularly oscillating Vc with large amplitude 
appears via cooperation of spikes in the fully-occupied stripes, and its global frequency /c 
(= 10.2 Hz) is smaller than that for the partially-occupied case (in fact, /c for the fully- 
occupied case is nearly equal to the most probable frequency of v for the isolated single 
case). 

To further understand the full and the partial occupations in the raster plots, we examine 
the local and the global output signals for both the fully-occupied and the partially-occupied 
cases of D = 20. Since our neural network is globally coupled, any local neuron may be 
a representative one. Unlike the isolated single case (shown in Fig. [T]l, each local neuron 
within the network is stimulated by a synaptic current which is directly related to the aver- 
age of firing events of all the other neurons. First, we consider the fully-occupied excitatory 
case, and investigate the time series of the local potential vi of the 1st neuron and the global 
potential Vc which are shown in Fig.[3la). The 1st neuron fires spikes phase-locked to Vc at 
every cycle of Vc (i-S-i the local potential vi exhibits a nearly regular phase-locked oscilla- 
tion). To confirm this 1:1 phase locking, we collect 5 x lO'* ISIs from all local neurons, and 
obtain the ISI histogram [see Fig.jBjb)]. A single peak with a small width is located around 
the average value (= 97.9 ms) which agrees well with the global period 7c of Vc- This 1:1 
phase locking is expected to occur for the fully-occupied case because the global frequency 
/c (= 10.2 Hz) of Vc is nearly equal to the most probable frequency (= 10.3 Hz) of noise- 
induced oscillation for the uncoupled single neuron. As a result of such 1 : 1 phase locking, 
full occupation occurs in the raster plot. Second, we study the partially-occupied inhibitory 
case. Figure [Sjc) shows the time series of the local potential vi of the 1st neuron and the 
global potential Vc. The 1st neuron exhibits intermittent spikings phase-locked to Vc at ran- 
dom multiples of the period Tc (=54.2 ms) of Vc. In addition to these intermittent spiking 
phases, hopping phases (exhibiting small subthreshold oscillations) appear in most of global 
cycles. We also note that after occurrence of a spiking, recovery from a hyperpolarized to 
a resting state is made during the next global cycle. Hence, a "preparatory" phase without 
spiking and hopping (for preparing for generation of the next spike or hopping) follows 
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Fig. 3 1:1 phase locking (spikings phase-locked to the global potential Vq at every cycle of Vq) in N{= 
10^) globally coupled excitatory ML neurons: (a) time series of the local potential V[ of the 1st neuron and 
the ensemble-averaged global potential Va and (b) interspike interval (ISI) histogram with a single peak. 
Stochastic phase locking leading to stochastic spike skipping (intermittent spikings phase-locked to Vg at 
random multiples of the period of Vq) in N{= 10^) globally coupled inhibitory subthreshold ML neurons: (c) 
time series of v'l and Vq, and (d) ISI histogram with multiple peaks. For both inhibitory and excitatory cases, 
Idc = 87 /iA/cm~, D = 20 jlA ■ ms'/^^/cm^, and J = 3 mS/cm^. Vertical dashed lines in (a) and (c) represent 
the times at which local minima of Vq occur. Each ISI histogram in (b) and (d) is made of 5 x lO"* ISls and 
the bin size for the histogram is 5 ms. Vertical dotted lines in (d) denote integer multiples of the global period 
Ta (=54.2 ms) of Vq- 



each spiking phase (see the gray parts). In this way, the local potential exhibits an irregular 
firing pattern consisting of randomly phase-locked spiking and hopping phases and prepara- 
tory phases. We note that a regular fast global oscillation with a small amplitude emerges 
from these irregular individual oscillations. To confirm the stochastic spike skipping (arising 
fromn stochastic phase locking) in the local potential, we collect 5 x 10^ ISIs from all local 
neurons and get the ISI histogram. Multiple peaks appear at multiples of the period Tc of 
the global potential Vg. However, due to appearance of preparatory cycles, the 1st peak of 
the histogram (which is the highest one) appears at 2 Tq (not Tc). Hence, local neurons fire 
mostly in alternate global cycles. Due to this stochastic spike skipping partial occupation 
occurs in the raster plot. Similar skipping phenomena of spikings (characterized with multi- 
peaked ISI histograms) were found in singl e noisy neuron models driven by a weak periodi c 
external force for the stochastic resonance ( lLongtin|[l995 l2000l : iGammaitoni et al. Ill998l) . 
Unlike this single case, stochastic spike skipping in networks of inhibitory subthreshold 
neurons is a collective effect because it occurs due to a driving by a coherent ensemble- 
averaged synaptic current. As discussed in details in Sec tion [H similar results on skipping 
were also obtained in simplified networks of IF neurons jBrunel and Hakirn ] | 19991 : iBmnell 
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Fig. 4 "Statistical-mechanical" spike-based coherence measure in N{= 10') globally coupled ML neurons 
for Idc = 87 /jA/cm^, D = 20 /iA ■ ms'/-/cm^, and 7 = 3 mS/cm^. Partially-occupied inhibitory case: (a) 
time series of the global potential Vg, (b) plot of the global phase (3>) versus time (<), and (c) plots of O, 
(occupation degree of spikes in the ;th stripe), (pacing degree of spikes in the ith stripe), and M, (spiking 
measure in the ith stripe) versus i (stripe). In (a) and (b), vertical dashed and solid lines represent the times 
at which local minima and maxima (denoted by open and solid circles) of Vq occur, respectively and G, 
(;' =1,2) denotes the ith global cycle. Fully-occupied excitatory case: (d) plots of O,, P,, andM, versus i. 



l2000l ; lBmner& Wang 11200 3V However, there are some differences in the spiking pattern of 
individual neurons and the effect of noise on the collective coherence. 

Our main purpose is to quantitatively characterize the stochastic spiking coherence 
which is well visualized in the raster plot of spikes. To measure the degree of stochas- 
tic spiking coherence seen in the raster plots in Fig. |2] for D = 20, we introduce a new 
"statistical-mechanical" spike-based coherence measure by considering the occupation 
pattern and the pacing pattern of spikes in the stripes of the raster plot. Particularly, the 
pacing degree between spikes is determined in a statistical-mechanical way by quantifying 
the average contribution of microscopic individual spikes to the macroscopic global poten- 
tial V(j. The spiking coherence measure M,- of the ;th stripe is defined by the product of the 
occupation degree (9, of spikes (representing the density of the ith stripe) and the pacing 
degree Pi of spikes (denoting the smearing of the ith stripe): 

Mi = Or Pi. (13) 
The occupation degree (3, in the zth stripe is given by the fraction of spiking neurons: 
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where A'; is the number of spiking neurons in the ith stripe. For the full occupation, O, = 1, 
while for the partial occupation O; < 1 . The pacing degree Pj of each microscopic spike in 
the ith stripe can be determined in a statistical-mechanical way by taking into consideration 
its contribution to the macroscopic global potential Vc- Figure |4l a) shows a time series of 
the global potential Vc; local maxima and minima are denoted by solid and open circles, re- 
spectively. Central maxima of Vc between neighboring left and right minima of Vc coincide 
with centers of stripes in the raster plot. The global cycle starting from the left minimum of 
Vc which appears first after the transient time (= lO-' ms) is regarded as the 1st one, which 
is denoted by Gi. The 2nd global cycle G2 begins from the next following right minimum 
of Gi, and so on. Then, we introduce an instantaneous global phase 4>(?) of Vc via linear 
interpolation in the two successive subregions forming a global cycle (Freund et al. 20Q^. 
The global phase 4>(f) between the left minimum (corresponding to the beginning point of 
the ith global cycle) and the central maximum is given by: 

4>(r)=2.(i-3/2) + .(^ ^J-^^^^^,.„^ j for r"<^<^f"'"' 0=1,2,3,...), (15) 

and 4>(f) between the central maximum and the right minimum (corresponding to the be- 
ginning point of the {i + l)th cycle) is given by: 

/ t — f^""-*' \ 

'f(0 = 2^(.--l) + W^-^^^ forf("'-^'<r<f(:f 0=1,2,3,...), (16) 

where is the beginning time of the ith global cycle (i.e., the time at which the left 
minimum of Vc appears in the ith global cycle) and is the time at which the maximum 
of Vc appears in the ith global cycle. The global phase <P in the first two global cycles 
is shown in Fig. |Hb). Then, the contribution of the kth microscopic spike in the ith stripe 

(s) 

occurring at the time to Vc is given by cos<Pk, where 4>j; is the global phase at the 

('>) 

kth spiking time [i.e., 4>ic = 4>(f^ )]. A microscopic spike makes the most constructive (in- 
phase) contribution to Vc when the corresponding global phase <Pt is 2Kn (h = 0, 1,2, . . .), 
while it makes the most destructive (anti-phase) contribution to Vc when 4>,' is 2K{n — 1/2). 
By averaging the contributions of all microscopic spikes in the ith stripe to Vc, we obtain 
the pacing degree of spikes in the rth stripe, 

1 ^' 

/^■ = - J^cos4>,, (17) 

k=i 

where Si is the total number of microscopic spikes in the ith stripe. By averaging M, of 
Eq. ( 113b over a sufficiently large number A'j of stripes, we obtain the 'statistical-mechanical" 
spike-based coherence measure M/. 

We follow 3 X lO-' stripes and get O,, P,, and Mj in each ;th stripe for the partially- 
occupied inhibitory case of Fig.|2jb). The results are shown in Fig. |4tc). For comparison, 
we also measure the degree of stochastic spiking coherence seen in the raster plot of Fig. [2d) 
for the fully-occupied excitatory case and give the results in Fig.|4ld). We note a distinct dif- 
ference in the average occupation degree (O,), where (• • •) denotes the average over stripes. 
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Fig. 5 Stochastic spiking colierence for various values of D in N{= lO') globally coupled inhibitoiy sub- 
threshold ML neurons when /oc = 87 flA/cm^ and 7 = 3 mS/cm^. (a) Raster plots of spikings (i: neuron in- 
dex) and plots of number of spikes (A'j) versus time (/) and (b) interspike interval (ISI) histograms for (al) and 
(bl) D = 10 /jA ■ ms'/Vcm^ (a2) and (b2) D = 20 /lA ■ ms'/Vcm^. (a3) and (b3) D = 30 ^A ■ ms'/'/cm^. 
and (a4) and (b4) D = 32 ^iA-ms'^^/cm^. The bin size for Ns in (a) is 1 ms. Each ISI histogram in (b) is 
made of 5 x lO'* ISls and the bin size for the histogram is 5 ms. Vertical dotted lines in (b) represent the 
integer multiples of the global period Tq of Vq; Tq = (bl) 67.4 ms, (b2) 54.2 ms, (b3) 48.6 ms, and (b4) 47.7 
ms. (c) Plot of (O,) (average occupation degree of spikes) versus logigZ). (d) Plot of (P,) (average pacing 
degree of spikes) versus logmZ). (e) Plot of M, ("statistical-mechanical" spike-based coherence measure) 
versus logjoD. To obtain (O,), {Pi), and M, in (c)-(e), we follow the 3 x 10^^ stripes for each D. Open circles 
in (c)-(e) denote the data for Z) = 20 flA ■ ms'/^/cm^. 



For the inhibitory case, partial occupation with very small (O,) (= 0.106) occurs due to 
stochastic spike skipping, while for the excitatory case full occupation with (O,) = 1 occurs 
as a result of 1:1 phase locking. For the partially-occupied case of inhibitory coupling, the 
average pacing degree (Pi) (= 0.766) is large in contrast to (O,), although it is smaller than 
(Pi) (= 0.91 1) for the fully-occupied case of excitatory coupling. As a result, the "statistical- 
mechanical" coherence measure (which represents the collective coherence seen in the 
whole raster plot) is 0.081 for the partially-occupied inhibitory case which is very low when 
compared with Ms (= 0.9 11) for the fully-occupied excitatory case. The main reason for the 
low degree of stochastic spiking coherence is mainly due to partial occupation. 

Finally, by varying D we characterize the stochastic spiking coherence in terms of the 
"statistical-mechanical" measure Ms in a population of inhibitory subthreshold ML neurons 
when loc = 87 and J = 3. Such stochastic spiking coherence may be well visualized in the 
raster plots of neural spikes. Figures [3a l)-|5ja4) show the raster plots and the temporal plots 
of the number of spikes (A'^,) for A' = 10^ in the coherent region for D = 10, 20, 30, and 
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32, respectively. The corresponding ISI histograms are also given in Figs.|5tbl)^5tb4). We 
measure the degree of stochastic spiking coherence in terms of (O,) (average occupation 
degree), (Pi) (average pacing degree), and Ms for 20 values of D in the coherent regime, and 
the results are shown in Figs. |5jc)^e). (Here, for comparison with other values of D we 
include the case of D = 20 which is studied in details above; open circles in Figs.[5lc)-[5je) 
represent the data for D = 20.) For the coherent case of D=20, the raster plot in Fig.|5ta2) 
consists of relatively clear partially-occupied stripes with (O,) = 0.106 and (Pi) = 0.766. 
Such partial occupation results from stochastic spike skipping of individual neurons seen 
well in the ISI histogram of Fig. |Hb2) with clear (well-separated) multiple peaks appear- 
ing at nTc {T(j'. global period of V(j and n=2,3,...). As the value of D is increased from 20, 
the average occupation degree (O,) increases only a little, as might be seen from the raster 
plots and the plots ofN^ in Figs.[5];a3) and|5];a4) ((O,) =0.114 and 0.115 for D = 30 and 
32, respectively). This slow increase in (O,) is well shown in Fig.|5jc). However,the aver- 
age pacing degree {Pi) for D > 20 decreases rapidly, as shown in Fig. |Hd). For example, 
stripes in the raster plots of Figs. [5ja3) and |3a4) become more and more smeared, and 
hence the average pacing degree is decreased with increase in D. This smearing of stripes 
can be understood from the change in the structure of the ISI histograms. As D is increased, 
the heights of peaks are decreased, but their widths are widened [see Figs.|5jb3) and|5jb4)]. 
Thus, peaks begin to merge. This merging of peaks results in the smearing of stripes. We 
note that for large D merged multiple peaks in the ISI histogram are directly associated with 
smeared partially-occupied stripes in the raster plot. Thus, for D > 20 the degree of stochas- 
tic spiking coherence is rapidly decreased as shown in Fig. [Sje), mainly due to the rapid 
decrease in {Pi) (a little increase in (O,) has negligible effect). Eventually, when passing the 
higher threshold DJJ (ci 33.4) stripes no longer exist due to complete smearing, and then 
incoherent states appear. 

By decreasing the value of D from 20, we also characterize the stochastic spiking co- 
herence. As an example see the raster plot of spikes in Fig.[5tal) for D = 10. The average 
occupation degree is much decreased to (O,) = 0.044, as can also be expected from the plot 
of Ns- This rapid decrease in (O,) for D < 20 can be seen in Fig.|5lc). Contrary to (O,) the 
average pacing degree {Pi) (= 0.684) for D = 10 is decreased a little when compared to 
the value of {Pi) (= 0.766) for D = 20; only a little more smearing occurs. For this case, 
the ISI histogram has multiple peaks without merging like the case of D=20, as shown in 
Fig l3bl). However, when compared to the case of /) = 20 the heights of peaks decrease, 
and the average ISI increases via appearance of long ISIs. Thus, the degree of stochastic 
spiking coherence for D = 10 (M, = 0.032) is much decreased mainly due to rapid decrease 
in the average occupation degree (O,) (small decrease in {Pi) has only a little effect). In fact, 
as can be seen in Fig.[5jd) there is no noticeable change in {Pi) for 10 < Z) < 20; as D is 
decreased from 20 to a value of D (:i 14), {Pi) increases slowly, and then it begins to de- 
crease. Thus, for D < 10 both (O,) and {Pi) decreases so rapidly, as shown in Figs.|5tc) and 
|5jd). Hence, as D is decreased from 10 the degree of stochastic spiking coherence decreases 
rapidly due to both effects of (O,) and {Pi). Eventually, when D is decreased through the 
lower threshold D*i (2± 9.4), completely scattered sparse spikes appear without forming any 
stripes in the raster plot, and thus incoherent states exist for D <D*i .In this way, we char- 
acterize the stochastic spiking coherence in terms of the "statistical-mechanical" measure 
Ms in the whole coherent region, and find that reflects the degree of collective coherence 
seen in the raster plot very well. When taking into consideration both the occupation degree 
and the pacing degree of spikes in the raster plot, a maximal spiking coherence occurs near 
£) ~ 20 [see Fig.|5te)]. As discussed in details in Section[T] Ms is in contrast to the conven- 
tional measures where any quantitative relation between (microscopic) individual spikes and 
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the (macro scopic) global potential Vg is not considered [e .g., the "thermodynamic" order 



parameter taolomb & Rinzel Pi 994; 'Hansel & Mato 
based measure (Wang & Buzsaki 1996; White et ah 
of precise spike timing dMainen & Seinowski 1119951 ; 



2003 ), the "microscopic" correlation- 
1998 ), and the PSTH-based measures 



Schreiber et aTlbOOBh l. 



4 Summary 



By changing the noise intensity, we have characterized stochastic spiking coherence in an in- 
hibitory population of subthreshold biological ML neurons. In a range of intermediate noise 
intensity, a regularly oscillating global potential Vc (with reduced amplitude and the in- 
creased frequency) emerges via cooperation of irregular individual firing activities. We note 
that this stochastic spiking coherence may be well visualized in the raster plot of spikes. For 
the case of a coherent state, partially-occupied stripes appear in the raster plot. This partial 
occupation occurs due to stochastic spike skipping, which is shown clearly in the multi- 
peaked ISI histogram. Recently, Brunei and Hakim (1999) have obtained similar results on 
spike skipping in a simplified network of integrate-and-fire (IF) neurons which are randomly 
connected via instantaneous inhibitory synapses (modeled by the delta function with a trans- 
mission del ay). [Such skipping was also reported in the network of inhibitory and excitatory 
IF neurons terunel Il200ol ; Brunei & Wang Il2003l) .1 For the case of the IF neuron a spike is 
triggered instantaneously when the membrane potential reaches a threshold, in contrast to 
the case of the biological ML neuron where dynamics of voltage-gated ion channels leads to 
spike generation. As the stimulus exceeds a threshold, the firing frequency of the IF neuron 
begins to increase from zero (i.e., the IF neuron exhibits the type-I excitability), unlike the 
case of the type-II ML neuron (used in our computational study). These type-I IF neurons are 
driven by random excitatory inputs activated by independent Poisson processes. In this sim- 
plified network, it was shown that weakly coherent global activities (roughly corresponding 
to the global potential Vg) emerge from irregular firing patterns of neurons [refer to Figure 3 
in the paper of Brunei and Hakim (1999)]. The collective coherence was well shown in the 
temporal auto-correlation function of the global activity. Although the basic results in the 
simplified network are similar to ours in a realistic network, there are some differences in the 
spiking pattern of individual neurons. The stochastic spike emission of individual IF neurons 
was well shown in the Poissonian ISI histogram showing little indication of the collective be- 
havior. Since no clear multiple peaks appear in the Poisson-like ISI histogram, the individual 
IF neurons exhibit more stochastic spike emission than the ML neurons, although weak col- 
lective coherence occurs in both cases. There are some additional differences in the effect of 
the external noise on the collective coherence. As the magnitude of the external noise is in- 
creased from zero and passes a threshold, the global activity w as found to exhibit a transition 
from an oscillatory to a stationary state llBmnel and Hakinii1 [r999). Thus, incoherent states 
appear for the case of strong noise. On the other hand, through competition between the 
constructive role (stimulating coherence between noise-induced spikes) and the destructive 
role (spoiling the collective coherence) of noise, stochastic spiking coherence (in our work) 
appears in a range of intermediate noise intensities. Hence, incoherent states appear for both 
cases of sufficiently strong and weak noise. These differences in both the spiking pattern 
of individual neurons and the noise effect on the collective coherence seem to arise mainly 
from different types of drivings on individual neurons; IF neurons are driven by random ex- 
citatory inputs activated by Poisson processes, while ML neurons are driven by a Gaussian 
white noise in addition to a subthreshold DC stimulus. The main purpose of our work is 
to quantitatively measure the degree of stochastic spiking coherence seen in the raster plot. 
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We have introduced a new type of spike-based coiierence measure by taking into con- 
sideration the occupation degree and the pacing degree of spikes in the stripes. Particularly, 
the pacing degree between spikes is determined in a statistical-mechanical way by quantify- 
ing the average contribution of (microscopic) individual spikes to the (macroscopic) global 
potential Vc- This "statistical-mechanical" measure is in contrast to the conventional 
measures (e.g., the "thermodynamic" order parameter, the "microscopic" correlation-based 
measure, and the PSTH-based measures of precise spike timing) where any quantitative re- 
lation between the microscopic individual spikes and the macroscopic global potential Vq 
is not considered. In terms of Mj, we have quantitatively characterized the stochastic spik- 
ing coherence, and foimd that reflects the degree of collective spiking coherence seen in 
the raster plot very well. Finally, we also expect that Mj may be implemented to character- 
ize the degree of collective spiking coherence in an experimentally-obtained raster plot of 
neural spikes. 
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